##############################################################################################################
# This script figures out, for each shape of type 1 (e.g. administrative area), and for each shape of type 2,
# the area of the intersection
#
# It assumes that both shape files are given in the same geometry and that this geometry is "flat" (e.g. Lambert93).
# It outputs a CSV file with that information.
# If the intersection is empty, no row is provided in the output file.
#
# Example use:
#	python3 getIntersectionAreas.py input1.shp shapeId1 input2.shp shapeId2 output.csv
##############################################################################################################

import sys, time, csv
import cartopy
import cartopy.io.shapereader as shpreader
from myutilities import printRemainingTime, appendData2File

def segmentsOverlap(seg1, seg2):
    # This function determines whether two segments overlap, meaning that they have a non-empty intersection.
    # seg1 should be a list or tuple giving: [min1, max1]
    # seg2 should be a list or tuple giving: [min2, max2]
    if seg1[0] >= seg2[0] and seg1[0] <= seg2[1]:
        return True
    if seg1[1] >= seg2[0] and seg1[1] <= seg2[1]:
        return True
    if seg1[0] <= seg2[0] and seg1[1] >= seg2[1]:
        return True
    return False

def boxesOverlap(box1, box2):
    # This function determines whether two boxes overlap, meaning that they have a non-empty intersection.
    # Here I refer to a box as an n-dimension extension of a rectangle:
    # 	in each of the n dimensions, there is a lower bound and an upper bound.
    # box1 should be a list/tuple of dimension n, where box1[ii] gives [min, max].
    # box2 should be similar.
    NumDims = len(box1)
    for dd in range(NumDims):
        if not segmentsOverlap(box1[dd], box2[dd]):
            return False
    return True

def getBoundsFromShapes(shapes, epsilon=0):
    # This function takes shapes as inputs and returns the bounds (minX1,maxX1,minX2,maxX2), so I can plot a map from them.
    # When calling shape.bounds, we get: [minX1, minX2, maxX1, maxX2]
    # This function returns the mins of the mins, and the max's of the max's (for both X1 and X2),
    # and possibly ands/removes an epsilon (for padding) around these bounds.
    Inf = float('inf')
    minX1 = Inf
    minX2 = Inf
    maxX1 = -Inf
    maxX2 = -Inf
    for shape in shapes:
        mybounds = shape.bounds
        minX1 = min(minX1, mybounds[0])
        minX2 = min(minX2, mybounds[1])
        maxX1 = max(maxX1, mybounds[2])
        maxX2 = max(maxX2, mybounds[3])
    minX1 = minX1 - epsilon
    minX2 = minX2 - epsilon
    maxX1 = maxX1 + epsilon
    maxX2 = maxX2 + epsilon
    return [minX1, minX2, maxX1, maxX2]

[inputShapeFile1, shapeId1, inputShapeFile2, shapeId2, outputFile] = sys.argv[1:6]

# Start output file
fout = open(outputFile, 'w')
fout.write(shapeId1 + ',' + shapeId2 + ',IntersectionArea\n')
fout.close()

# Read administrative shapes
shapes1 = list(shpreader.Reader(inputShapeFile1).records())
shapes2 = list(shpreader.Reader(inputShapeFile2).records())


# 1) Calculate the bounds of all admin shapes
bounds1 = []
for s1 in shapes1:
    bounds1.append(getBoundsFromShapes([s1.geometry]))

# 2) Calculate the bounds of all coverage shapes
bounds2 = []
for s2 in shapes2:
    bounds2.append(getBoundsFromShapes([s2.geometry]))

# 3) Iterate
box1 = [[], []]
box2 = [[], []]
output = []
NumShapes1 = len(shapes1)
NumShapes2 = len(shapes2)
tic = time.time()
for idx1 in range(NumShapes1):
    # First check whether the rectangles defined by their bounds overlap
    s1 = shapes1[idx1]
    box1[0] = [bounds1[idx1][0], bounds1[idx1][2]]
    box1[1] = [bounds1[idx1][1], bounds1[idx1][3]]
    for idx2 in range(NumShapes2):
        s2 = shapes2[idx2]
        box2[0] = [bounds2[idx2][0], bounds2[idx2][2]]
        box2[1] = [bounds2[idx2][1], bounds2[idx2][3]]
        # If so, calculate area. If not zero, append it to results
        if boxesOverlap(box1,box2):
            shapesIntersection = s1.geometry.intersection(s2.geometry)
            myarea = shapesIntersection.area
            if myarea > 0:
                output.append([s1.attributes[shapeId1], s2.attributes[shapeId2], myarea])
    if (idx1+1) % 20 == 0:
        printRemainingTime(tic, idx1+1, NumShapes1)
    if (idx1+1) % 10 == 0 or idx1 == (NumShapes1-1):
        appendData2File(output, outputFile, ',')
        output = []
